Granular Brownian motion 



A. Sarracino, D. Villamaina, G. Costantini and A. Puglisi 

CNR (ISC-SMC) - p.lc A. Moro 2, 00185, Roma, Italy 

Dipartimento di Fisica, Universita Sapienza, p.le A. Moro 2, 00185, Roma, Italy 
E-mail: ale . sarracino@gmail . com,dario . villamainaQromal . infn. it 

Abstract. We study the stochastic motion of an intruder in a dilute driven granular 
gas. All particles are coupled to a thermostat, representing the external energy source, 
which is the sum of random forces and a viscous drag. The dynamics of the intruder, 
in the large mass limit, is well described by a linear Langevin equation, combining 
the effects of the external bath and of the "granular bath" . The drag and diffusion 
coefficients are calculated under few assumptions, whose validity is well verified in 
numerical simulations. We also discuss the non-equilibrium properties of the intruder 
dynamics, as well as the corrections due to finite packing fraction or finite intruder 
mass. 
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1. Introduction 

Granular materials in the fluidized state [U [2] have represented, during the last 10-15 
years, an excellent benchmark for new and old theories in non-equilibrium statistical 
mechanics: the presence of non-conservative forces make unavailable the standard 
tools used at equilibrium, such as Gibbs measure, equipartition, thermodynamic limit, 
Einstein relation and more [3j HJ E]. It is therefore necessary to resort to more 
fundamental theories, from the Boltzmann equation up to stochastic processes and 
modern generalizations of statistical mechanics to non-equilibrium states [TJ, [8] . 

In order to achieve a stationary state, the fundamental ingredient is an external 
source of energy, required to compensate the energy lost in inelastic collisions. The 
role of energy source is played by some injection mechanisms, depending upon the 
experimental setup, e.g.: a box with a vibrating wall, a layer (or more than one) placed 
over a vibrating plate, a gas flux going through orifices in the box walls, etc. The 
different mechanisms may produce quite different states with different symmetries: for 
instance a layer over a vibrating plate is homogeneous on average, while a boundary 
driving (e.g. a shaken box wall) leads to spatial gradients and currents j9j [101 ITT] . 

From the point of view of a tracer particle, however, the dynamics is always 
of a similar kind: the tracer interacts, in a random sequence, with the surrounding 
particles and with the energy source. The ratio between frequencies of interaction 
dictates the relevance of tracer-particle collisions with respect to exchanges between 
the tracer and the source. Of course, in a boundary-driven setup, the statistics of 
collisions suffered by the tracer depends upon the distance from the energy source. 
Anyway, the random motion performed by the granular tracer should always take into 
account the two contributions: collisions with other granular particles and interaction 
with the energy source. In more idealized setups, the so-called Homogeneous or non- 
Homogeneous Cooling States, no energy injection is involved: in this regimes, anyway, 
a collisional stationary state cannot be achieved and experimental verification is very 
difficult to be achieved. 

Here we consider a model commonly used in the theoretical literature on granular 
fluids: all grains are coupled to a thermostat-like energy source, with a typical interaction 
time Tb which is usually taken larger than the inter-particles collision time r c [T2], [13] - The 
stationary granular gas obtained in this way, is then used as a "granular bath" where 
a massive intruder performs a non-equilibrium Brownian motion, still being coupled 
to external energy source. The result is a double bath whose properties are analyzed 
starting from a linear Boltzmann-Lorentz-Fokker-Planck equation, which is treated in 
the diffusional approximation (large mass) to be cast into a Langevin equation. 

Self-diffusion of an intruder [TJJ [15] or a tracer [161 EE] has been previously studied 
in the Homogeneous Cooling State. The same calculations have been performed for 
models with an impact velocity dependent restitution coefficient [18] and a good review 
of main results can be found in the textbook [8] . Diffusion in a stationary granular fluid 
obtained by imposing shear boundary conditions has also been considered [19] . 
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Up to our knowledge this is the first time that self-diffusion of a large mass intruder 
is studied for a model with homogeneous energy injection, considering explicitly the 
effect of a "double bath" , i.e. of both sources of noise, granular and external respectively. 
A Langevin equation (I3T]) with expressions for the tracer temperature, Eq. ( 134)) . mobility, 
Eq. ( 1321) . and diffusion coefficients, Eq. (1351) . all involve the interplay of both energy 
sources. The large mass limit, together with the Molecular Chaos assumption (due to 
diluteness) guarantees that a granular Fluctuation-Dissipation theorem holds, where 
the ratio between diffusion and mobility is simply given by the intruder granular 
temperature [201 121] • 

In Section [2] we introduce the model (granular gas, thermostat and intruder); in 
Section[3]the Kramers- Moyal expansion and the large mass limit are discussed, leading to 
the Langevin formulation in Section 13.21 Numerical experiments (Molecular Dynamics 
and Direct Simulation Monte Carlo) are performed to study the limits of the used 
assumptions in Section [4] and finally conclusions and perspectives are drawn in Section [5j 



2. The model 



We consider a gas of N granular spheres in d dimensions, each sphere has index i, with 
i G [1 , AT] , and mass m,. Particle 2 = 1 (referred to as "the intruder") has mass M 
and radius R, while all other particles (usually denoted as "the gas" ) have mass m and 
radius r. The parameter e = ^Jm/M will be used for large mass expansion. The system 
is contained in a box of volume V = L d , much greater than the volume occupied by the 
particles, so that the hypothesis of molecular chaos applies. We denote by n = N/V the 
density of the gas and by <fi the occupied volume fraction (in d = 2 it is, for instance, 
<J) = 7t[(N -l)r 2 + R 2 ]/V). 

The intruder and the gas particles undergo binary instantaneous inelastic collisions 
when coming at contact, with the following rule 

v, = w[ (i + a ) [U - v') • a] a (1) 



where Vj (v^) and (v^) are the post and pre-collisional velocities of particle % (particle 
j), respectively; a G [0,1] is the restitution coefficient 0, and a is the unit vector 
joining the centers of the colliding particles. The mean free path of the intruder 
is Iq = l/(n(r + i?) d_1 ). Two kinetic temperatures can be introduced for the two 
species: the gas granular temperature T g = m{vf)/d (i > 1) and the intruder granular 
temperature T tr = M{v 2 )/d. 

In order to maintain a fluidized granular gas, an external energy source is coupled 
to every particle in the form of a thermal bath. The motion of a particle i with velocity 
Vj is then described by the following stochastic equation 

mm(t) = -76 V, (t) + fi(t) + £ b (t). (3) 

| for simplicity we consider the restitution coefficient to be equal for all particles. 
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Here fj(i) is the force taking into account the collisions with other particles, and £ b (t) is 
a white noise, with (£&(£)) = and (£&,ia(^)£&,j,8(^ / )) = ^T b ^ h 5ij5 a p5{t — t'), where Latin 
indices refer to particle labels while Greek indices denote Cartesian coordinates 0. 

The effect of the external energy source balances the energy lost in the collisions 
and a stationary state is attained. Several temporal scales are important in this system: 

• r| , the mean free time between collisions of a gas particle; 

• r* r , the mean free time between collisions of the intruder; 

• r b 9 = m/'jb the typical interaction time of the bath with gas particles; 

• r b tr = M/^fb the typical interaction time of the bath with the intruder. 

When 76 is small enough to have the mean free times t§ and r* r smaller than the 
interaction times r b and T b r , inelasticity is sufficient to put the gas out of equilibrium: 
this is reflected, among other things, in the failure of equipartition T g < T b and T tr < T b . 
It is also known that T g ^ T tr (221 123] . 

The main goal of this note is to show that, in the limit of large mass M, the 
force fi acting on the intruder can be expressed by means of a Langevin-like formula 
fi(£) = — 7 9 V(£) +£ g (t), providing explicit expressions for 7 3 and (£ g (i)£ g (t')). 

In order to do that, let us start by writing the coupled Boltzmann equations for 
the probability distributions P(V,t) and p(v, t), denoting (for simplicity) with V and 
v the intruder velocity and the gas velocity, respectively 
dP(V,t) 



dt 
dp(v, t) 



J dV'[W tr (V\V')P(V',t) - W tr (V'\V)P(V,t)} + B tr P(V,t) 

m f dv'[W g (v\V) P (V,t)-W g (V\v)p(v,t)]+B gP (v,t) 

+ J[v\p,p], ( 4 ) 
where Bt r and B g are two operators taking into account the interactions with the thermal 
bath. In these equations the effects of the collisions for the tracer and the gas particles 
are described by, respectively, 

Wir(V|V) = x J dV J dap(V, t)Q [- (V - v') ■ a] (V - v') • a 



6® l 



e 2 



V - V + -^—^{l + a) [(V -V)-a]a\ (5) 



and 



W g (v\V) = §JdV J daP(Y, t)Q [- (V - V) • a] (V - V) ■ a 

x tf« |v - v' + — + a) [(V - V) • a] a j , (6) 



where Q(x) is the Heaviside step function, S^(x) is the Dirac delta function in d 
dimensions, and x = 92 ^ R ^ > Qzi? + R) being the pair correlation function for a gas 

§ We use a constant 7&, but in principle this coefficient may depend on the mass and on the radius of 
the particle, since it is only a model description of a more complicate interaction with plates, walls or 
fluids going through the granular medium. 
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particle and an intruder at contact; in the expressions §5§ and (jB]) we have assumed 
that the probability P 2 (|x — X| = r + R,*V , v, t) that a collision between the intruder 
and a gas particle occurs, when they have velocities V and v and positions X and x 
respectively, is given by the Enskog approximation [8] 

P 2 (|x - X| = r + R, V, v, t) = g 2 (r + R)P(V, t)p(v, t) (7) 

which is a small correction to Molecular Chaos, taking into account density correlations 
near the intruder; the terms describing the action of the thermal bath read 

B tr P(V, t) = ^ A [V P(V, t)} + ^A F [P(V, t)] (8) 

Mv,*) = -|" [vp(v,t)] + ^A„[p(v,t)], (9) 
mov m 

where is the Laplacian operator with respect to the velocity; finally, the Boltzmann 
collision operator for the particle-particle interactions J[v|p, p], can be found in many 
papers, see for instance [23]. In view of the fact that it is not relevant for the rest of the 
paper, we omit its explicit expression. 

2.1. Decoupling the gas from the tracer 

The two Boltzmann equations appearing in the system (j3J) are coupled through the 
terms involving Wt r and W g . Nevertheless, if the number N of granular particles is large 
enough, the term W g can be neglected because of the factor 1/N in Eq. ([6]). Hence, the 
surrounding gas is weakly perturbed by the tracer and fast and homogeneous relaxation 
is expected. One assumes that the probability distribution function p(v) is stationary 
and, following numerical evidence (verified below) it is approximated with a Gaussian 
function with variance T g /m: 



p(v) = — = exp 

y/(2nT g /mY 



mv 2 



2T g j 



(10) 



Substituting Eq. ( flOl) into Eq. (jSJ), and projecting the velocities along the collision 
direction and the orthogonal one, the integral can be solved [25], yielding 

1 



W tr (V'\V)=xk(e)- 2 (V;-V c ^ 



y/2irT g /m 

x exp{-m[fc(e)- 1 (K-^) + K] 2 /(2T 5 )} 



11^ 



where V a = V-a (note that a is parallel to V — V) and k(e) = (1 + a)e 2 /(l + e 2 ). From 
now on we specialize to the two dimensional case, where the above equation simplifies 
to 

1 



W tr (V'\V) = X 



^2nT g /mk(e) 2 

x exp {-m[VZ -V a + k(e)V a } 2 /(2T 3 A;(e) 2 )} . (12) 
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As discussed in details below, once the gas is decoupled from the intruder, the dynamics 
of the tracer alone is Markovian, and it is known that such transition rates satisfy 
detailed balance with respect to a Gaussian invariant probability -P(V) [25] (the 
temperature of the tracer, in that case, where m = M, is given by fz^T g [26]). 

2.2. Granular temperature of the gas 

The granular temperature T g can be obtained from the Langevin equation Q. Indeed, 
multiplying by v(t) and averaging, one gets 

^m|(v 2 (t)) = - 7 ,(v(t) 2 ) + (v(t)f(t)) + <v(t)£ 6 (f)>. (13) 

At stationarity, the l.h.s. of the above equation vanishes and (v(£)£ 6 (£)) = I^T^jm. 
The term (v(t)f(t)) represents the average power dissipated by collisions, which we 
assume to be dominated (this is true for N large enough) by gas-gas collisions: 

(v(t)f(t)) = -{AE) coh (14) 

where AE = l/8m(l — a 2 )[(v! — v 2 ) • a] 2 is the energy dissipated per particle and the 
collision average is defined by 

(■■■)coi = Xgj da J dvi J d\ 2 ■ • .p(vi)p(v 2 )6[-(vi - v 2 ) ■ ct]|(vi - v 2 ) ■ a\ 

where \g = and l 9 = l/(n(2r) d ~ l ) is the mean free path for gas-gas collisions and 
g' 2 (2r) is the pair correlation function for two gas particles at contact. The integral in 
Eq. (fT4|) can be computed by standard methods [8], and, in two dimensions, yields 

(AE U = X3 ^pf>T^. (15) 

Substituting this result into Eq. (fl3~|) and recalling that T g = m(v 2 )/2, one finally obtains 
the implicit equation 



27fe 

which can be solved to obtain T g . 

3. Kramers-Moyal expansion for the tracer-gas collision operator 

With the assumption discussed above, the system of equations d3J is decoupled. That 
allows us to write the following linear Master Equation for the tracer 

dP{ ^ t] = L gas [P(V,t)] + L bath [P(V,t)}, (17) 

where L gas [P(V, t)] is a linear operator which can be expressed by means of the Kramers- 
Moyal expansion [27] 

L 9as [P(V,t)] = £ -1 j ° Dg,JV)P(V,t), (18) 
n=l Ji " " " J« 
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(the sum over repeated indices is meant) with 

"/"I,. (V) = ^ / dV'(l£ - ^) • • • (VI - V jn )W tr (V'\V), (19) 

and W tr is given by relation ( 1T2"|) . The second term in the Master Equation represents 
the interaction with thermal bath: 

L bath [P(V,t)] =B tr P(V,t). (20) 

In the limit of large mass M, i.e. small e, we expect that the interaction between the 
granular gas and the tracer can be described by means of an effective Langevin equation. 
In this case, we keep only the first two terms of the expansion 



L gas [P(V,t)} = -^[A (1) (V)P(V,t)] + t*[D®{V)P(y,t)]. (21) 

A justification of this truncation, in the limit of small e, comes from observing that 
terms ■ are of order e 2n : this can be obtained by plugging Eqs. ([T|) (for the case of 
the tracer, i.e. V = Vi) into (Tl9~l) . 



It is useful at this point to introduce the velocity-dependent collision rate and the 
total collision frequency 



r(V) = J dV'W tr (V'\V), (22) 
to = [ dV P(V)r(V). (23) 



The former quantity can be exactly calculated, giving 

— /rp \ 1/2 



X 



(eg + 2)/ ( — I + e g ^ I — 



(24) 



where the rescaled variable q = V/^T g /M is introduced in Appendix through 
Eqs. (1A.14P and I n (x) are the modified Bessel functions. To have an approximation 



of oj, on the other side, one has to make a position about P(V). Let us take it to be 
a Gaussian with variance T tr /M. The consistency of this choice will be verified in the 
following section. With this assumption, the collision rate turns out to be 



u = X V2^^T g /m + T tr /M = X V^ (j^j ' W 1 + ^e 2 = u K(e), (25) 



T 
1 <j 



1 /2 

where u = xV^ and K(e) = ^Jl + ^e 2 . 

3.1. Large mass limit 



We are then able to compute the terms Z? t - and appearing in L gas . The result 
and the details of the computation of these coefficients as functions of e are given 
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in Appendix. Here, in order to be consistent with the approximation in ([21]), from 
Eqs. (IA.15j) we report only terms up to 0(e 4 ) 

= - x ^\ x {\ + a )e 3 + 0(e 5 ) 
m 

T \ 1/2 

J- n 



i\ { i + a y Vx + (e 5 ) 

1 J 

= -u (l + a)e 2 V x + O(e 5 ) (26) 
£>« = - w (l + a)e 2 V y + 0(e 5 ) (27) 

(rp \ 3/2 
(l + a)V + 0(e 5 ) 

= ^(l+«)¥|0( f 5 ) (28) 

D$ = 0(e e ). (29) 

The linear dependence of Dp upon Vg (for each component (3), allows a granular 
viscosity 

r] g = w (l + a)e 2 . (30) 

In the elastic limit a — > 1, one retrieves the classical results: 7/ g — >■ 2woe 2 and 
£)£z = D^/ — > 2uj e 2 ^. In this limit the Fluctuation-Dissipation relation of the second 
kind is satisfied [281 EHj , i.e. the ratio between the noise amplitude and j g , associated 
to the same source (collision with gas particles), is exactly T g /M. When the collisions 
are inelastic, a < 1, one sees two main effects: 1) the time scale associated to the drag 
r g = 1/rjg is modified by a factor ijp, i.e. it is weakly influenced by inelasticity; 2) the 
Fluctuation-Dissipation relation of the second kind is violated by the same factor ^jp. 
This is only a partial conclusion, which has to be re-considered in the context of the full 
dynamics, including the external bath: this is discussed in the next section. 

3.2. Langevin equation for the tracer 

Putting together the results in Eqs. ( 12611291) with Eqs. ( TlTll2~TT) . we are finally able to 
write the Langevin equation for the tracer 

MY = -TV + 8, (31) 

where r = 7& + 7 9 and 8 = £ b + £ , with 

7 3 = Mr/g = Mu (l + a)e 2 = w„(l + a)m (32) 

'1 + a r 
~2 

concluding that the stationary velocity distribution of the intruder is Gaussian with 
temperature 

r t = gllM. (34) 

Ib + lg 



lbT b + 7 9 — — T g 



SM-t'), (33) 
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Equation ( 131 j) is consistent with the Gaussian ansatz used in computing oj . Note that 
the above expression for T tr is consistent with the large mass expansion obtained in 
Eqs. ( 1291) only if it is dominated by T g , for instance when 7 9 ^> 75 (see discussion at 
the end of Appendix A ). In the opposite limit, the tracer dynamics is dominated by 
the coupling with the external bath and the typical velocity of the tracer cannot be 
taken sufficiently small with respect to the typical velocity of gas particles, making the 
expansion unreliable. In this case, however, if the diameter of the intruder is similar to 
that of the gas particles, it is reasonable to expect similar collision frequencies: the gas 
particles will therefore be dominated by the external bath and the whole system will be 
very near to equilibrium [30, EG] • 

For the self-diffusion coefficient it is immediately obtained 

D tr = [ dt(V x (t)V x (0)) = ^ = ^^ff^ - (35) 

Solving numerically the equation (fTB"j) and substituting the result into the above 
equation, one can study D tr as a function of the restitution coefficient a (this is done 
numerically in the next section). When all other parameters are kept constant and a is 
reduced from 1, the behavior of D tr is non-monotonic, it decreases, has a minimum and 
then increases for lower values of a. Anyway, this minimum is expected for quite low 
values of a or high values of the packing fraction 0, where the approximations involved 
in this theory are not good. For this reason, at the values of parameters chosen to have 
a good comparison with simulations, this non-monotonic behavior is not observed. 

It should be also noticed that, in the Homogeneous Cooling State, the self-diffusion 
coefficient at a given granular temperature increases as a is reduced from 1, i.e. it has 
an opposite behavior with respect to the present case [Ml [15]. Other studies on different 
models of driven granular gases have found expressions very close to Eq. ( )32|) . which is 
not surprising considering the universality of the main ingredient for this quantity, i.e. 
the collision integral [3 [31] . 



3.3. Energy fluxes and detailed balance 

A few comments are in order, at this point, concerning the non-equilibrium properties of 
this system. The first question comes about the term ^-jp which multiplies T g in Eq. (134|) . 
It is easily explained with the following argument [25]: we have assumed that the tracer 
feels no memory of past collisions, which means that any post-collisional correlation 
with recoiling gas particles is lost. With these assumption, the fate of recoiling particles 
can be ignored and the dynamics concerns only the intruder: 

V = V>-(l + a) J ^{(V>-v)-d-]a, (36) 

where v is the pre-collisional velocity of the colliding gas particle (randomly extracted 
from the given distribution p(v)). Then, one simply observes that for any value of a, 
M and m, such rule can be rewritten as an elastic collision rule with an effective mass 
M' = 2 A f~ Hm - mK — iW for large intruder mass. This is equivalent to say that the 
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tracer has elastic interactions with the gas particles, with an effective mass M', and 
therefore feels an effective temperature of the gas T' g = jpT g = ^MpT 9 . Note that this 
argument, for m = M, gives the formula T' = ^^T g , which has been derived for the 
first time in [26] . 

The energy injection rates of the two thermostats [32] are 

Q b = <V(f) • - 76 V)> = 2^(T b - T tr ) (37) 
Q 9 = (V(t) • (£ 9 - 7g V)> = 2^(T' g - T^) (38) 

It is easy to see that the balance of fluxes Qb = —Q g is equivalent to formula ( j34"|) for 
T tr - This balance implies that, if T tr < T, then T tr > T' . When a < 1, the two fluxes 
are different from zero, i.e. energy is flowing from the external driving, through the 
tracer, into the granular bath. 

Apparently, this contradicts the "equilibrium" nature of the Langevin equation ( I3TT) : 
the tracer dynamics is Markovian and stationary, and the equation satisfies detailed 
balance with respect to the Gaussian invariant distribution. As already discussed 
in [25], this is not a paradox but only a consequence of Molecular Chaos and the 
decoupling assumption which allows us to write Equation (TIT]) : here we have employed 
the Enskog approximation, which is a weak modification of Molecular Chaos, still 
preserving Markovianity, i.e. no memory terms appear in Eq. (J3J). The absence of 
memory implies that both and £ 5 are white noises and makes them indistinguishable: 
an observer which can only measure V(i) cannot obtain separate measures of Qb and 
Q g , but only a measure of the total energy flow Q = M(V ■ V) = which hides out 
the presence of energy currents. A more detailed analysis, e.g. by relaxing the Enskog 
approximation, should put in evidence the different time-correlations of the two baths: 
eventually, the observer, by means of some "filter", should be able to sort out their 
different contributions Qb and Q g . This is an interesting example where memory plays 
a crucial role in the non-equilibrium characterization of a system [33] . 

We expect that time reversibility (detailed balance) is a symmetry, for the intruder, 
which is broken in the following cases: 1) at small values of M (this is different from the 
case discussed in [25], where the intruder was not in contact with the external bath); 
2) when the non-Gaussian behavior of the gas velocities is taken into account; 3) when 
the tracer has asymmetric properties with respect to some spatial axis [M]; 4) when 
Molecular Chaos (or its weak Enskog correction) is violated [5]. 

4. Numerical simulations 

In this Section we report the results of Molecular Dynamics (MD) simulations of 
the model, together with Direct Simulation Monte Carlo (DSMC) simulations [35] 
incorporating the Enskog correction, and compare them with our theoretical predictions. 
In all simulations we have kept constant the dimension d = 2, the mass of gas particles 
m = 1 and the radii r = R = 0.01, as well as the properties of the bath TJ, = 1 and 
7b = 0.1; instead we have varied N, M, a and <j) (values of L and n can be obtained 



Granular Brownian motion 



11 



from the knowledge of r and 0). We have used the Carnahan-Starling expression for g 2 
at contact [36]: g 2 {r + R) — (1 — — < / ) ) 2 - For the chosen values of < 0.07, it 

is always g 2 {r + R) < 1.12. In all simulations we have also checked that the Gaussian 
approximations for the velocity distributions of gas particles and for the intruder are 
satisfied, observing very small values for the second Sonine coefficient a 2 < 0.02 [8]. 

In Fig. [T]we show the velocity-velocity autocorrelation function C(t) = (V x (t)V x (0)) 
of the tracer for different values of its mass M = 100, 25, 5, 2 in a dilute and moderately 
inelastic case: a = 0.8, and = 0.00785 (and N = 10 4 for MD). We can clearly observe 
that in the case of large mass M = 100 the Langevin equation (131!) describes very well 
the dynamics of the tracer. Indeed, in that case, the numerical results are consistent 
with the theoretical prediction 

C(t) = ^e-iK (39) 

As expected, for smaller values of M, the numerical results move away from the 
analytical ones and large corrections to the exponential decay do appear. The deviations 
are observed (and are quantitatively similar) for both MD and DSMC results, implying 
that they are due, as expected, to the breakdown of the large mass expansion, rather 
than that of Molecular Chaos. For MD results we have noticed that, going from N = 10 3 
to N = 10 4 , the comparison with DSMC (and with theory at large M) is improved. 

In order to check the validity of the hypothesis of molecular chaos, we report the 
results of MD and DSMC simulations for higher packing fractions in Fig. [2j keeping 
M = 100, N = 10 4 (in MD) and a = 0.8: since the clean part of the decay of C(t) 
is always exponential, we focus only on the two parameters of interest, i.e. T tr and j g . 
One clearly observes that, increasing the packing fraction, the discrepancy between the 
theoretical value and the values obtained from MD, increases. On the other side, DSMC 
always gives results very close to theory, as expected. The Enskog approximation (J7|), 
which does not take into account memory effects, is no longer valid in MD at high packing 
fraction, while always holds in DSMC. In order to enforce this statement, we computed 
the following correlation coefficient: Cy u = (&v x Sum) where we introduced the 

stochastic variable u m {t) given by the averaged x-component velocity of the particles 
lying, at time t, in a fixed area around the tracer. In particular, 8V X and 5u m measure 
the deviations of V x and u m from the average values, which tend to for a large number 
of measures. The coefficient defined above must be zero, if molecular chaos holds; on 
the contrary we observed that its value sensibly increases as the packing fraction gets 
higher. For example for = 0.00785 ,Cy Um = 0.005 whereas, for = 0.07, Cy Um = 0.07. 

Finally let us compare the diffusion coefficient D tr = J °° dt C (t) measured in MD 
and DSMC with the theoretical value obtained through Eqs. (|T6|) and (!35l) . In Fig. [3] 
we show our results at different values of a, keeping fixed M = 100, = 0.00785 and 
iV = 10 4 (in MD). Again there is a perfect match for DSMC, while MD simulations 
present a small discrepancy which becomes more evident at small values of a. We have 
again verified that this discrepancy is a finite effect and is reduced as N increases. 
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Figure 1. (Color online). The autocorrelation function C(t) = (V x (t)V x (0)) is 
measured in MD and DSMC simulations (black circles and red diamonds, respectively) 
for M = 100, 25, 5, 2 in the model with restitution coefficient a = 0.8 and packing 
fraction <fi = 0.00785 and coupled to a thermal bath with 7b = 0.1 and Tf, = 1. The 
blue lines show the theoretical predictions of Eq. (|39l) . 



5. Conclusions 

While many papers have been devoted to the large-mass diffusive properties of an 
intruder in a cooling granular gas, the driven case, somehow, has received less 
attention [19]: this is in contrast with the fact that, in real experiments, the most 
common situation is a driven granular gas. The problem, at the level of the basic 
assumptions treated here (Enskog approximation, negligible non-Gaussianity and large 
separation of timescales between collisions and driving), does not pose particular 
conceptual difficulties, nevertheless it reveals to be already quite rich. The external 
driving mechanism, characterized by a temperature T& and the "internal" granular bath 
at temperature T g < T b , sum up together in giving a linear Langevin dynamics for 
the intruder, provided that the collision frequency between the intruder and the gas 
particles is larger than the frequency of interaction with the bath. Such Langevin 
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Figure 2. (Color online). The temperature Ttr (top panel) and the drag coefficient 7 9 
(bottom panel) measured in MD (black circles) and DSMC (red diamonds) is plotted 
for different values of the packing fraction <j> = 0.00785, 0.2, 0.5, 0.7 in the model with 
AI = 100, N = 10 4 (in MD) and a = 0.8 (error bars fall within the symbols). The 
dashed blue lines show the theoretical predictions following from Eqs. (|34l32p . 




0.6 0.7 0.8 0.9 



a 

Figure 3. (Color online). The diffusion coefficient of the tracer D tr is measured in 
MD (black circles) and DSMC (red diamonds) simulations for different values of the 
restitution coefficient a = 0.9, 0.8, 0.7, 0.6 in the model with M = 100 and packing 
fraction </> = 0.00785. The dashed blue line shows the theoretical prediction following 
from Eqs. (|16I35[) . In the inset the same curve is plotted in the whole range a € [0, 1]. 



equation predicts for the "intruder temperature" T tr a weighted sum (with weights 
given by the drag coefficients of the two baths) of TJ, and T' g = ^-^-T g , i.e. the intruder 
feels the surrounding gas to be at a different temperature T' g < T g , because of non- 
conservative interactions. The self-diffusion coefficient is even more interesting, showing 
a non-trivial non-monotonic behavior with a minimum at low values of the restitution 
coefficient. Our results lose validity when the mass of the intruder is reduced, when the 
packing fraction of the gas is increased, when the inelasticity is too low to disregard 
non-Gaussian corrections, and when the interaction times of the two baths become 
comparable. 
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It is interesting to discuss what is happening at moderately high packing fractions 
4> ~ 10%: we have seen that the Enskog approximation is not very good to predict the 
intruder dynamics, because it is missing memory effects mediated by the surrounding 
fluid. A scenario which can be conjectured is the following: the gas may display 
two typical relaxation times, a local one related to collisions r re \ ~ r| and a global 
one r' Tel > r re i, which is due to diffusion of slower modes (e.g. hydrodynamics). If 
Kei > T c' > T reh one nas that the intruder feels a "locally equilibrated" surrounding 
granular gas. In this case it is reasonable to replace Eq. ( jTUl) with 



where u and T g are some local velocity and temperature fields which change on timescales 
larger than r* r (and correspondingly large spatial scales). A partial verification of this 
scenario has been mentioned at the end of [5], but requires further investigation. 
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Appendix A. Calculation of first two coefficients of the Kramers-Moyal 
expansion 

For larger generality (whose motivation is discussed in the Conclusions), in this 
Appendix we discuss the case where the gas surrounding the intruder may have a non- 
zero average u[JJ: 




(40) 





which is a simple task involving only the definition of new shifted variables 



c = 



c = 



V-u 



V'-u. 



(A.2) 
(A.3) 



We are interested in computing 




x exp j-m [c' a + (fc(e) - 1) C(7 ] 2 /{2T g k(e) 2 )} 



(A.4) 



|| note that in all the cases discussed in the main text, we have always taken u = 0. 
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Figure Al. An example for the change of variables (c' x ,c' y ) — > (c CT , c^.), introduced in 
Eq. (|A.5[) . Such change of variable, when inverted, has two possible determinations: 
in this example both represented vectors c' yield the same (c CT ,c^). 



In order to perform the integral, we make the following change of variables (see 
Fig. IA1I for an example) 



c' — c T 



C(T c x 



d = c' 



C y c y 



C T . c x )^ ~\~ (c,. Cy 



C-r C-x) 2 ~\~ (c„ Cy 



"y ^y) 



c y - c v 



c x - c x y + {d y - Cy y v ^{d x - Cx y + {d y - Cy y 



which implies 



dc' = dc' x dc' —> dc a dc' a \J\, 



where 



\J\= A ^ • r 



■ 2 -c 2 ) 



is the Jacobian of the transformation. The collision rate is then 

^2 

2 



r(V) = x 



7T 



2T g /m 



[c 2 + 2T g /m)I 



mc 



+ c 2 h 



mc 



4T„ 



where I„(x) are the modified Bessel functions. For D,- we can write 



/ + OO /*00 1 

dcj dc'M-c,)\J\-== 
-oo </<:„ x/Z-Kln 



/mk{e)' 



x exp {-m [c' a + (fc(e) - l)c ff f /(2T 3 A;(e) 2 )} 



+c 



rfc^ / dc'^c'i - Ci 



c„ - c n 



(A.5) 

(A.6) 
(A.7) 

(A.8) 



V c2 - c l 

^--exp {-m [c' a + (k(e) - l)c a } 2 / (2T g k(e) 2 )}A.9) 



'2i:T g /mk{e) 
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where we have enforced the constraint of the theta function, namely c a G (— c, +c), with 
c = \/ c x + c y- Notice that the integral in dd a is lower bounded by the condition c' a > c a 
which follows from the definition of c a . In order to compute the integral, we have to 
invert the transformation ( 1A.5|) . That yields two determinations for the variables c' x and 
c' y (see Fig. ETJ) 



(A) 



c -^- (c a c x + c y Sign(c x )^c 2 - c 2 a 
c x Sign{c x )^Jc 2 - c 2 a 



C y c y 



c a c x - c y Sign(c x )^c 2 - c 2 a 
c a c y + c x Sign(c x )^c 2 - c 



Then the integral flA.9[) can be written as 



1 i*c /»oo 

^ 1} (V) = - / dcj dd a [(<£ - C x )^ + (4 - C x )^] \J\ 
L J-c Jc a 

exp \-m [ C ; + (*(e) - l)c a f /(2T g k(e) 2 ) 



y/2irT g /mk(t 



(A.10) 



yielding 



3lr 



mix 

2Tn 



' 4T„ 



(2) 

Analogously, for the coefficients D\- } one obtains 



(c 2 + 3T 9 /m)/ (^) + (c 2 + T g /mM^) 



-9 
J. 



(c 2 + 3T g /m)I (^) + (c 2 + T a /m)hC^r) 



(A.H) 



1 1 



2/ 



J-c Jc a L J 

1 



X 



w '2TTT g /mk(e) 2 
1 1 A;(e) 2 /2mvr 

2T ^TV~rr e 



exp -m [<£ + - l)cj' /(27>(e) 2 ) 



c 2 (4c 2 + c 2 ) + 3T 5 (7c 2 + 3c 2 )/m + 15T 2 / 



^ 2 'm 2 



mc 



4T„ 



+ 



-c 2 + c 2 



c 2 (4c 2 + c 2 ) + T 9 (13c 2 + 7c 2 y )/m + 3T 2 /m 2 % " 



mc 



4T n 



(A.12) 
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C x~ C 



y ^y> 



ic> x -c x r\c> y -c y )^ 



]\J\ 



X 



^2TiT g /mk(e) 2 
1 1 k(e) 2 



exp [-m [c' a + (Me) - l)c CT ] 2 /(2T,A;(e) 2 )} 



2 in 



T 



9 Ct.Ch 



(c 2 + 4T 9 /m)/ 



Then we introduce the rescaled variables 



mc 



+ 



c 4 + 2c 2 T g /m - 2T 2 /m 



mc 

4i; 



(A.13) 



m 



in 



(A.14) 



obtaining 



^(V)= -\^-\ x k{e)ee-^ 
3 Iq V 2 m 



2^2 



( e V + 3)J (if )+(eV + l)/ 1 



2„2 



3/2 



[e 4 g 2 (4g 2 + q 2 y ) + 3e 2 (7g 2 + 3q 2 y ) + 15] J 



x 



+ 



qI + qI 



eV(4g 2 + g, 2 ) + e i (13g 2 + 7g 2 ) + 3 



1 1 1 



T 



3/2 



'eV + 2e 2 g 2 -2 



(eV + 4)/ 



4 



e 2 g 2 



h 



e 2 q 2 
4 



(A.15) 



Up to this last results we have not introduced any small e approximation. The next 
step consists in assuming that q ~ 0(1) with respect to e, which is equivalent to assume 
that c 2 ~ Tg/M: this assumption must be compared to its consequences, in particular 
to Eq. (1341) ; the assumption is good for not too small values of a and for 7 9 3> 7&, i.e. 
when T tr ~ T g . When this is the case, expanding in e and using that Iq{x) ~ 1 + x 2 /A 
and Ii(x) ~ x/2 for small x, one finds Eqs. (T29]) . 
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